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Abstract 

An outstanding issue in the computational analysis of time dependent prob- 
lems is the imposition of appropriate radiation boundary conditions at artificial 
boundaries. In this work we develop accurate conditions based on the asymptotic 
analysis of wave propagation over long ranges. Employing the method of steepest 
descents, we identify dominant wave groups and consider simple approximations to 
the dispersion relation in order to derive local boundary operators. The existence 
of a small number of dominant wave groups may be expected for systems with 
dissipation. Estimates of the error as a function of domain size are derived under 
general hypotheses, leading to convergence results. Some practical aspects of the 
numerical construction of the asymptotic boundary operators are also discussed. 


1 Introduction 

Many interesting and important problems involving wave propagation are posed on un- 
bounded spatial domains. Examples include wave propagation in the atmosphere and 
the ocean as well as flows past bodies and in channels and tubes. For purposes of numer- 
ical computation, an artificial boundary is often introduced. For long time computations 
the interaction of the solution and the artificial boundary cannot be avoided. This may 

‘Supported, in part, by Space Act Agreement C-99066-G and by NSF Grant No. DMS-8905314. 


1 



result in unacceptably large errors throughout the computational domain, especially if 
the system has instabilities. (For an example of this see Hagstrom and Keller [11].) 

The primary purpose of this work is to develop accurate boundary conditions to be 
imposed at such boundaries. We are also interested in establishing error estimates and 
convergence theorems, standard components of theoretical numerical analysis which seem 
rarely to have been studied in the context of time dependent partial differential equations 
on unbounded domains. For linear, separable differential operators, the exact boundary 
conditions may be represented in terms of appropriate eigenfunction expansions and 
transform variables. (See, e.g., Gustafsson and Kreiss [7].) For example, suppose x is 
the spatial coordinate normal to the boundary and the x-dependence of the transform 
solutions takes the form e Xl ^ x with s the dual variable to time and / indexing a tangential 
normal mode. An exact relation at the boundary is given in transform space by: 



ui — 0 . 


(1) 


Unfortunately, the expression of this relation in the original variables is typically nonlocal 
in both space and time. For computational efficiency, the added storage and arithmetic 
operations required by the implementation of nonlocal conditions must, to the extent 
possible, be avoided. Ideally, local boundary operators would be used. These may be 
obtained using polynomial or rational approximations to the dispersion relation, A/(s), 
which in turn can generally be accurate only in a restricted neighborhood of transform 
Space. 

An approach to the derivation of boundary conditions is, evidently, to identify re- 
gion^) in transform space where polynomial or rational approximations are to be made 
and, then, to compute the coefficients of the approximation. In their pioneering study 
of hyperbolic problems, Engquist and Majda [5] considered a particular high frequency 
limit. For problems with dissipation, on the other hand, this limit is less likely to lead 
to accurate results. In this work we consider the use of asymptotic expansions of waves 
propagating over long distances computed using the method of steepest descent. A 
consequence of the dissipative terms in the equations studied here is the association of 
growth or decay with each wave group. By locating minimum decay (maximum growth) 
rates, we locate appropriate regions for the required approximations. That is, we iden- 
tify a small number of dominant wave groups, characterized by (5,-, /j), and compute local 
linear approximations to the dispersion relation, 


A(,(s) sa A; t (S,) + Aj.(S,)(s — s,). 


(2) 


Substituting this into (1) leads to a local operator. The complete asymptotic boundary 
condition is defined by the composition of a small number of these operators: 





Si) U = o. 


( 3 ) 
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In Sections 2 and 3, we present in detail the asymptotic analysis and the subsequent 
derivation of asymptotic boundary conditions. Some numerical considerations are dis- 
cussed in Section 4. In Section 5 we derive estimates of the error as a function of the 
size of the computational domain. These lead to the convergence of the solution of the 
problem on the truncated domain to the solution on the full domain. 

A primary motivation of this work is the development of accurate boundary condi- 
tions at artificial boundaries for the Navier-Stokes equations. This is carried out in [10] 
for incompressible flows, where extensive numerical experiments are described. Earlier 
applications of some of the ideas given here appear in [8] and [9]. Boundary conditions 
for similar equations with constant coefficients have been derived by Halpern [12] and 
Halpern and Schatzman [13]. 


2 Asymptotic Expansions 


We consider, for definiteness, a general system of equations in a semi-infinite channel: 


du du du d 2 u „d 2 u 

~at +u Tx + v Ty + Wu = A dT> +B w' 


( 4 ) 


a-' > 0, yo<y<y\- (5) 

These are supplemented by boundary and initial conditions defining a signalling problem: 

u{x,y, 0) = 0, (6) 

D 0 u(x,y 0 ,t) = 0, (7) 

D\u{x,yx,t) = 0, (8) 

E o u(0,y,t) = g(y,t). (9) 

We assume that the matrices Do , D\ and Eq are such that the problem is well-posed. 
We allow a stratified medium; that is U, V , W, A and B are functions of y. Note that 
equation (4) may be a far field approximation to a problem whose coefficients are either 
nonlinear or functions of x. 

A representation of the solution of problem (4-9) may be obtained by means of Laplace 
transforms and eigenfunction expansions. The eigenvalue problem to be solved is: 


svi + A iUvi -f V + Wv t = XfAvi + 2/o < i/ < 2/i, 

D 0 vi{y 0 ]s) = 0, D 1 vi(y l ;s) = 0. 


( 10 ) 

( 11 ) 
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For solutions, it, which grow at most exponentially in time we may restrict attention to 
eigenvalues, A/, satisfying: 

3fc(A i(s)) < 0, 3fc(s) sufficiently large. (12) 

We denote by M the set of indices of eigenvalues which meet the condition above and 
will refer to the function A i(s) as the dispersion relation. Let 

J roo . . 

' e~ st 9(yJ) dt - ( 13 ) 

0 

The assumption of well-posedness then implies the existence of a unique collection of 
functions c/(s) such that: 

Eo[j2ci{s)v,(y,s)\=g{y,s). (14) 

V&V ) 

If ci(t ) is the inverse transform of cj(s), a final expression for u may be obtained: 

u(x,y,<) = ^2ui(x,y,t), (15) 

1 6 AT 

t>t 

u,{x,y,t)= C[(p)qi(x,y,t — p)dp, (16) 

Jo 

where 

qi(x, y,t) = J c e’ <+Ai(s)x ^(y; s ) ds ( 17 ) 

and C is an appropriate inversion contour. 

To compute asymptotic expansions of u, valid for x large, we must evidently find 
expansions of q,. If (17) is evaluated along rays t = 'yx , x > 1, the exponent becomes 

x (7 s + A/(s)). (18) 

In order to use the method of steepest descent, we seek points s* such that: 

AJ(0 = -7, ( 19 ) 

»( 7 ) > 0 , ( 20 ) 

$(7) = 0. (21) 

Then, assuming that for 0 < 7 min < 7 < 7 max < 00 there exists $*( 7 ) satisfying (19-21) 
with inversion contour, (7, which can be deformed to the steepest descent path, we have. 

„(*,«,<) ~ e~(*'(-)<«* +A| ('-(*))) (22) 
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( 23 ) 


7 min ^ ^ 7 max- 

X 

Substituting these into (16) formally yields an approximation of u t for t > 7minS; 

ui(x,y,t)~ [ ci(p)<f)i(x,y,t — p) dp. (24) 

This representation has a simple interpretation: the signal data, c/, generates wave pack- 
ets which propagate at their group velocity. At the point (x,t), x > 1, the solution is 
approximately the superposition of waves generated at times varying from t - ^ max x for 

the slowest waves to t - 7 m i„x for the fastest. 

Now consider the specialization of these results to hyperbolic systems. In particular, 

we take the wave equation written as a first order system: 

a = b = w = 0 , i/=(J J), r=(o _°i). < 25 > 


and suppose that y 0 = 0 and y x = 1. For an appropriate choice of boundary conditions 

we have: 

A,(s) = - y/s ' 2 + / 2 7r 2 , (26) 

dXi —s (27) 

ds \/s 2 + / 2 tt 2 

From (27) we see that group velocities ranging from 0 to 1 are associated with values of 
s = *w, \u\ > 1 7T. Furthermore, »(A) = 0, that is the wave packets do not decay expo- 
nentially as they propagate. These observations hold in general for the high frequencies 

of all hyperbolic problems. (See, e.g. [15].) 

For problems with dissipation, on the other hand, it may be possible to further sim- 
plify the results. Then, some exponential decay rate may be associated with each wave 
group. That is: 

+ A ,(«•)) 4 0. (28) 


For general signal data the large x behavior will be dominated by the wave group with 
least decay (which may be growth for problems with instabilities). Therefore we seek 7 
such that the expression above is maximized. Setting to zero the derivative of the decay 
rate with respect to 7 yields: 


»(.• + ^(7 + A ',)) = 0 , 


(29) 


which by (19) reduces to: 


s(s‘M) = o. 


(30) 
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That is, assuming (19) defines a curve in s space, critical points of the decay rate occur 
as the curve crosses the imaginary axis. 

Suppose, for simplicity, that a unique solution of (30), 7/, exists. (If s* is imaginary, 
its complex conjugate must also be used.) So long as this wave packet is excited by 
the initial data, we expect that the dominant contribution to the convolution integral 
defining ut(x,y,t) will come from a neighborhood of 


(t~P) - 

-^- = 71. 

Introducing a local approximation to <j>i we obtain: 


(31) 


ui(x,y,t) ~ / 

J rr 


t 1mxr\X t ~ T> 


(32) 


Here we have: 


F(z, t) = e 4A,r 


[— 47 T/ 


— 47rA^r 


^ = lis* + A/(s*), 

d A, 

ds 


\ 1 d\^ . 

A i ~ ~r( s ) = ~ 7 /, 


(33) 


(34) 

(35) 



Id 2 A, 
2 ds 2 


(O, 


(36) 


where s* is evaluated at 7,. Also we assume that t > 7 t x . The expressions above 
lepresent a restriction to the neighborhood of a single point in the dispersion relation. 
In what follows this restriction will enable us to find an asymptotic boundary condition 
which consists of local operators. Furthermore, techniques will be given for the numerical 
computation of the various quantities defined in (34-36). 


3 Construction of the boundary conditions 


We now suppose that an artificial boundary is located at a: = r. The representation of 
the solution (15-16) may be manipulated to yield a variety of exact relationships at the 
boundary. For example, if the matrix A is positive definite, a characterization of the 
exact boundary conditions is that a collection of functions, r t (t), exists such that: 


£( t >V>0 

u(r,y,t) 


) = Y, / r ‘(p) ( 

/ leM Jo V 


I p ) 
qi(T,y,t - p) 


dp. 


(37) 
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In many cases the unknown functions r; may be eliminated to yield a direct relationship 
between u and J^. Of course, this condition will be nonlocal in y and t and, in general, 
too difficult to use. 

If, however, the asymptotic expansion given in (32) is valid, it may be used to develop 
a local asymptotic boundary condition. The Laplace transforms of and qi are related 

by: as 

■^(x,j/;s) = A ,(s)q,(x,y,s). 

The steepest descent result involves the restriction of the transforms to a neighborhood 
of $*(77). An asymptotic expansion of the x derivative may then be obtained by replacing 
Xi(s) by its Taylor series about the critical value of s: 

A i(s) « A,(s*( 7,)) + A }(s - s'( 7,)) + A *(s - s'( 7/)) 2 + ■ ■ . . (39) 



Using, for example, the first two terms we have: 


dqi 

dx 






(40) 


These may be substituted into (37) to finally obtain a condition on u /. The time derivative 
is brought outside the integral to further simplify the expression. This involves the 
neglect of terms from the limits of integration which should be exponentially small. The 
asymptotic boundary condition we propose is, then, given by: 


du t 

dx 



(41) 


A hierarchy of conditions may be obtained by use of more terms in the Taylor series. 
These would involve derivatives of higher order. Their stable implementation would 
require the use of Pade approximants, as discussed by Engquist and Majda [5]. For 
example, a potential approximation at the next order is: 



The condition derived involves only one normal mode, u/. We may, however, apply it 
directly to u if / is such that ^(Aj 1 ) is maximized. For problems where a small number of 
modes have similar minimum decay rates, a product boundary condition is used. Also, 
if the critical s m is imaginary we must include its complex conjugate in the product. We 
have in general: 




(complex conjugate) 


u = 0 . 


(43) 


We have never used more than two modes. 
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4 Numerical considerations 


In order to apply the asymptotic boundary conditions derived above, we must find 7/ 
as well as the first two terms of the Taylor series of Aj. For fixed s , the substitution 
w = Xv in conjunction with an appropriate discretization of the y derivatives leads to a 
generalized matrix eigenvalue problem: 

M(s)r = \Lr, (44) 


This may be solved using standard linear algebra software which implements, say, the 
QZ algorithm [6]. An analysis of the discretization error is given by Kreiss [17]. The 
asymptotic expansion requires that s be imaginary and ^ be real. As s varies along the 
imaginary axis, the latter implies a maximum or minimum for 5i(A). Our problem is, 
then, to maximize 3?(A/(s)) as s varies over the imaginaries. This is a line search problem 
for which many strategies have been developed, though they may be expensive to carry 
out. (See, e.g., Dennis and Schnabel [3].) For many of our examples, s — 0 has been the 
solution. A necessary condition for this is that jj(0) be real and negative. This may be 
checked by solving (44) only once. 

Once the critical value of s has been found, we must compute A/ . Differentiating (44) 
with respect to s yields: 




dM 

ds 


)r. 


(45) 


Since the matrix on the left is singular, A* = ^ may be obtained from a compatibility 
condition. If r is a left null vector we have: 


A 


1 

1 


fTdM r 

as 

r T Lr 


(46) 


The cost of this computation is typically negligible in comparison with the cost of finding 
the critical frequency. Expressions for more terms in the Taylor series can be similarly 
found. 

There are a variety of reasonable numerical implementations of the boundary con- 
ditions. For conditions involving derivatives of at most first order, many stable dis- 
cretizations are known. Product boundary conditions (for wave equations) are studied 
by Higdon [16]. He develops the useful principle of employing products of stable dis- 
cretizations. We have successfully employed this procedure throughout our numerical 
experiments. 
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5 Error estimates 


Estimating the error caused by the introduction of an artificial boundary can be broken 
up into two parts. The first is to estimate the residual resulting from the application of the 
boundary condition to the exact solution. The complete error estimate then follows from 
estimates of the solution of the initial boundary value problem in terms of inhomogeneous 
boundary data. The latter is simply a proof of well- posed ness. 

For reference we write down the problem satisfied by the error, e(x,y,t). 


Problem 1 

Le = 0, (x,y) e (0,r) x (?/o,J/i), (47) 

e(x,y,0) = 0, (48) 

s= 0, i = 0,1, (49) 

E o e(0,y,t) = 0, (50) 

Be(r,y,t) = Bu(r,y,t). (51) 


Here, L is the differential operator appearing in (4) and B is the asymptotic boundary 
operator appearing in (43). 

A general approach to the investigation of the effect of the boundary conditions on the 
well-posedness of an initial-boundary value problem is to freeze coefficients at each point 
of the boundary and to study solutions of the frozen coefficient system. In particular, we 
must show that no solutions of the frozen coefficient problem of the form: 

e K<+^+,*ay^ > 0, 3fty > o, (52) 

are in the null space of the boundary operator. For problems with decay in the dominant 
wave groups, that is < 0 in (43), it is clear that an eigensolution satisfying the 
conditions above cannot exist as: 


Aj - Aj/c) > 0. (53) 

(Recall that A/, is real and negative.) 

Translating this condition into a bonafide proof of well-posedness requires further 
assumptions on the coefficients of (4). Indeed, the conditions we have constructed may 
result in an overdetermined problem if A is rank deficient. Complete discussions are given 
by Eidel’man [4] for the parabolic case and by Strikwerda [19] for incompletely parabolic 
systems, though in general their results require boundary operators of low order. See 
also the more stringent requirement of dissipativity introduced by Barry, Bielak and 
MacCamy [1], which may be necessary for the estimates used below. We will simply 
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assume that Problem 1 is well-posed. We introduce appropriate Sobolev norms in the 
interior, || • ||(o,r)x(yo,»i). and on the boundary, || • || (#0 , Vl) and define: 


and, for functions of t: 


We then make: 



|w(*^ •, 0||(O t T)x(yo,yi)<ft, 


jttf 


jX — ||u^(x, *, f)||(y 0t J / i)df, 

roo 

Mi = / MOM*- 

Jo 


(54) 

(55) 

(56) 


Assumption 1 B C > 0 independent of r such that, 


Mil < C\\Bu\U, T . 


(57) 


(In what follows we mean all constants to be independent of r unless otherwise stated.) 

The requirement that C be independent of r can easily be dropped, though we must 
be able to estimate its growth in order to estimate the rate of convergence of the solution 
on the truncated domain. For example, in the error analysis of asymptotic boundary 
conditions for second order scalar hyperbolic equations in exterior domains as given by 
Bayliss and Turkel [2] and Hariharan and Hagstrom [14], algebraic growth of C with r is 
encountered. For the problems under consideration, however, Assumption (1) is typically 
satisfied due to the exponential decay of solutions. An interesting possibility is to use 
weighted norms in x, in which case C may decrease with increasing r. 

We proceed to estimate ||5u||i )T . We make a variety of simplifying assumptions, as 
our main purpose is to extract the dependence of the error on r as well as to give the flavor 
of the necessary computations. Many of the assumptions could be relaxed somewhat, 
though their verification for any practical problem is likely to be difficult. We Ireely quote 
error estimates from the asymptotic theory of integrals with large parameters. These may 
be found in any of a variety of texts on the subject, for example Sirovich [18]. 

We make, then, the following assumptions: 

Assumption 2 For all l € Af and all 7 G [0, 00) there exists a unique (up to complex 
conjugation) s*( 7) satisfying (19) and the contour of integration in (17) can be deformed 
to the steepest descent path through s * without passing through singularities of the inte- 
grand. Furthermore, for r sufficiently large, D*(qi — <f>i) is absolutely integrable in 7 and 
the integral is uniformly bounded in l where D represents a derivative in t or x and j < 2. 

Assumption 3 There exists /o (E Af such that a unique 7/ 0 exists at which 3?(7S* + 
A i 0 (s m )) = po attams a global maximum. If s*( 7; 0 ) = 0, then, for l = l 0 , ^(^*(7)) = 0 for 
all 7 in a neighborhood of 7/ 0 . Furthennore, 36 > 0 such that 9?(7<s* + A/(s’)) < Pq — 8 
for all 7 and l ^ Iq. 
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Assumption 4 There exists z such that, for 0 < j ^ 2. 

V [°° e z ^ p) hij(p)dp < I<0 < OO, ( 58 ) 

I# Jo 


pj(p) = 5R(5*(p)p + Aj(s (p))), 

= |»|(-.V(P))W ( |,.( P)| ,- + |A,(s»)|'). 

\/W(*-)l 

wi// always take r > z.) 

Assumption 5 There exists K\ > 0, I< 2 > 0 and, I< 3 > 0 such that, 

|c/|i < K\ |c; 0 |x, l / *0, 

||uj 0 ||l,T > -^2| c Joll|kjolll. T ’ 

|c, 0 |i < ifaIMko- 


(59) 

(60) 


(61) 

(62) 

(63) 


Assumptions (2 - 4) are constraints on the coefficients of equation (4). They are easily 
verified in the interesting special case of the advection diffusion equation with constant 
coefficients and the usual (e.g. Dirichlet or Neumann) conditions on the channel wall. 
Assumption (5) is in effect a restriction to signal data which excites the dominan wave 
group. If the data is such that only a certain wave group is excited, that information 

should be used to modify the asymptotic analysis. 

The main estimates are derived in the following collection of lemmas, whose proofs 

depend on the validity of Assumptions (2 - 5). 


Lemma 1 i 

\\Bqi 0 \\l,r = 0{~) • Ikolll.T- ( 

Proof: We have, by our assumptions on the properties of the transforms, 

Mi, = II*, Ml + °<7))- (65) 

Writing out the integral expression for IP 1 J 1 , and introducing the change of variables 
p — - yields: 

f e T f ^\ cos (rg(p) + 0(p))\h(p)dp. (66) 

V 2it Jo 

f(p) = ^t(s*(p)p+^i 0 ( s (?)))> 
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( 68 ) 

(69) 


g(p) = 3(s'(p)p + AJi-fp))), 
k(p) = ibkLi (p))lliw.wii 

vfwi ’ 


e{p) = --arg(A"(s*(p))). 


(70) 


*r K r SU T ’ } ' as a maximum at p = 7/0 = Po- The approximation of 

s integral by Laplace s method is slightly complicated by the presence of the absolute 

value of the oscillatory term. There are two cases to consider. If s*(p 0 ) = 0 then by 

Assumption (3) it is 0 in a neighborhood of p 0 and by direct computation we find that 

9 }fl r_“ ^tII ?u aCe ’ S f °T Ula may be dirCCtly appIied ‘ If 5 *(Po) ^ 0 then 
9 tf>oj - ^{s ) ^ 0. Then the generalization of Laplace’s formula given in Theorem 12) 

in the appendix applies. In each case the result is: 


II A Ikr = Ke T ">(l + o(-l)), 

T 2 


(71) 


where I< is independent of r. Similarly, asymptotic expansions of B( 
computed. In particular we have: 


dx 


hi o can be 


\\Bq lo \\ h T = \\B</> lo \\ hr (l + 0(-)). 

T 

As above we consider the integral expression for ||5^/ 0 (|: 

£ l er/<P) l c °s(t£(p) + 6(p))\h{ p )dp. 

Here / and g are as above while 0 and h are determined by: 


(72) 

(73) 

(74) 

Hp) = arg(5(s*(p),A, 0 (s*(p)))) - | arg (A£( s *(p))). ( 75 ) 

Again we will compute an asymptotic approximation to this integral using Laplace‘s 
formula if s*(p 0 ) = 0 or Theorem (2) if it isn’t. We note that 5 has been choL to have 
a simple zero at p — p 0? which implies that: 


Mp) = IK (■; »»)|| ( „, „, 


l g (C(p),A<„(a-()>)))| 


h(p) = 0(p- p 0 ), p y po # 

We thereby conclude: 

IlSfclkr = 0(i)|p 4 ||,, r . 

The conclusion of the lemma follows directly. 


(76) 

(77) 
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Lemma 2 


ll^tl Jo || llT =0(-)||« lo || 1 ,r. 


Proof: By direct computation, 

Bu lo = c; 0 * Bqt Q + boundary terms. (79) 

The asymptotic analysis indicates that the boundary terms are exponentially small. We 
then have, using Lemma (1) and Assumption (5), 

|| Biii 0 1| l lT « \\c l0 * Bq l0 \\ur 

5: |cfc,|i||^«olli.T (80') 

< ©(iJkMkfelll.r v ; 

< o(i)iiuj| 1>T . 

Lemma 3 There exists 7 / > 0 such that: 

P(« - = Ofe-niKII,.,. (81) 

Proof: We directly estimate ||£?(u — w; 0 )||: 

||B(u - « 4 )|„ » || E » * < A'IqJ, E WBMU.r, (82) 

Wo Wo 

where we have used Assumptions (2) and (5). From the integral representation of ||<^/|| 
we h ave ^ 

E l|B*||.,r < Jf E C eTMP> ~ h M d P< ( 83 > 

Wo V Z7r Wo J0 

where 

h(p ) = &(s*(p)p + M s *(p)))> ( 84 ) 

h,( P ) = M-;« W)lh-) N /| A ,, (s . (f>)) | — (85) 

However, by Assumptions (4) and (3), the expression on the right is bounded by: 




We have, therefore, for some constant K 4 : 




\\B(u-u lo )\\ hT <^e- ST (e T ^\c lo U). 
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By Assumption(5) and the asymptotic expansion of \\q to || computed in the proof of Lemma 
(1), the term in parentheses is bounded by ||u/ 0 || 1(T . We have thus shown that 7 ] > 0 can 
be chosen such that: 

||B( U - U ,„)||=0(e-’')K|| 1 , T , (88) 

completing the proof of the lemma. 

Lemma 4 There exists K such that: 

e Tfi ° 

ll-S^lll.r A — — ||u||i i0 . (89) 

Proof: We have, by lemmas (2) and (3): 

II < ||Bu,Ji, r + ||£(u- U ,„)||i, <0(i)IKIk,- (90) 

T 

Furthermore, using the asymptotic expansion of ||<^ 0 || 1>T computed in the proof of lemma 
(1) along with Assumption (5) we find: 

IKo||i, t < /^|Qoli||*olli,r < /i6e TMo ||u||i )0 . (91) 

Combining these equations yields the statement of the lemma. 

Finally, by combining Assumption (1) with lemma (4) we obtain the desired error 
estimate: 

Theorem 1 Suppose e is the solution of Problem 1 and Assumptions (1-5) hold. Then 
for r sufficiently large 3 K independent of r and the data such that: 

e T uo 

Nil < K— ||u||i, 0 . (92) 


6 Concluding Remarks 

We have developed a general technique for the construction of asymptotic boundary 
conditions for problems with dissipative wave propagation. Theorem 1 clearly implies 
the convergence of the solution on the truncated domain as r — ► oo if fj, 0 < 0. The 
factor of which is a direct consequence of the use of our asymptotic boundary con- 
ditions, is of importance when fio is small. Precisely such a situation occurs when our 
technique is applied to the incompressible Navier-Stokes equations for moderate to large 
Reynolds numbers. This is studied in detail in [10]. Applications to other problems in 
computational mechanics are also under consideration. 
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A Asymptotic expansion of an integral 

In this section we compute the asymptotic expansion as r — ► oo of the following integral 
which is needed for the derivation of our error estimates: 

1— I e T ^ p *| cos (rp(p) + 9(p))\h(p)dp. (93) 

Jp- 

The assumptions we will make are: 

Assumption 6 • The function f attains a global maximum at p Q G (p_,p+) and is 

thrice continuously differentiable in a neighborhood of p 0 with f'(po) ^ 0. 

• The function g is twice continuously differentiable in a neighborhood of po and 

g'(Po) ~f~ 0- 

• The function h is bounded outside a neighborhood of po and near p = po satisfies: 

h(p) = h 0 (p - p 0 ) m + MP " Po) m+1 (1 + 0(p - po)) , (94) 

where ho ^ 0 and, if m is odd: 

+ (95) 

• The function 9 is continuous in a neighborhood of po- 

• // p± = ±oo then f(p) — ► — oo at least algebraically as |p| — » oo. 
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The smoothness and decay conditions can be relaxed somewhat. If m is odd and (95) 
does not hold we must simply consider higher order terms in the expansion of / and h 
Despite its apparent simplicity, this integral does not seem to be discussed in the 
standard references. Due to the presence of the oscillatory term, Laplace’s method cannot 
be directly applied. Intuitively, we expect to obtain the dominant behavior by replacing 
| cos (rg + 0)| by its mean value. This is established in the following theorem. Although 
we have been able to estimate the contribution of the oscillatory terms, we have been 
unable to compute it to leading order for general / and g. 


Theorem 2 If Assumption (6) holds then as t — > oo: 




(96) 

where, 

f m even 

tm — | rrt odd ’ 

(97) 


r h 0 , m even 

Am “ | (hi - (m + 2)ho{$$), rn odd 

(98) 

Proof: We begin 

by replacing | cos (rg + 0)| by its Fourier series. 



2 00 

| cos 9\ = - + Cn cos 2n0 > 

* n=l 

(99) 

4 (— l) n+1 
~ n 7T 4n 2 — 1 

As the series converges uniformly, I may be expressed as the sum of integrals: 

(100) 


v 0 00 

I = / + -e T ^ p) h(p)dp + 5Zc n / n , 

Jp- K n= 1 

(101) 

where 

I n = f + e T/(p) cos (2 n(rg(p) + 0(p)))h(p)dp. 

dp- 

(102) 


The asymptotic expansion of the first integral follows directly from the use of Laplace’s 
method and yields the results stated in the theorem [18]. What remains is to show that 
the contribution of the oscillatory integrals is of lower order. Therefore we consider the 


asymptotic analysis of / n . 
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In the usual way we restrict the integral to a small, fixed interval about p 0 , introducing 
an error which is exponentially small. Using the fact that g'{p 0 ) ^ 0 we are able to make 
a change of variables so that the remaining integral becomes: 

CT/(P0) L C ° S (103) 

where 

F(u) = F 2 u 2 + F 3 u 3 + o(u 3 ), F 2 > 0, (104) 

and 

H{u,t) = H{u){l + 0(^)), (105) 

H(u) = H 0 u m + H lU m+1 + o(u m+1 ), H o ^0. (106) 

We break this integral into two parts, 

In = + /’), (107) 

II = J e~ Tp ' 2u2 cos (nTu)H 0 u m du, (108) 

ll = J e~ rF2u2 cos (titu)A(u, r)du. (109) 

Here we have, 

A(u,r) = ( e ~ T{F(u) ~ FiU ^ _ 1 )H(u,t) + ( H{u,t ) - H 0 u m ). (110) 

The asymptotic analysis of /* may be carried out using the method of steepest descent. 
The exponent has a single critical point at it = ~r leading to an exponentially small 
contribution. To estimate ll we first rewrite A: 

A(u,r) = ru TO+3 £ 1 (u, r) + ti m+1 <5 2 («, r), (HI) 

where 6 l and S 2 are bounded functions. Introducing the change of variables, 

v = ^/tF 2 u , ( 112 ) 

we obtain: 

9 m-f 2 f V + _ 2 T 

/ » = T 2 l e V cos (Jjrnv)(v m+ 3 d l (v,T) + v m+ 1 d 2 (v,T))dv. (113) 

Ileie, d y and d 2 are uniformly bounded and v± = F\Jt F 2 t. For convenience we extend 
the domain of integration to the entire real line, extending di and d 2 so that they remain 
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uniformly bounded. As the integrand is absolutely integrable on the entire real line 
uniformly in t, the additional term is o(r 2 ). To estimate the remaining integral we 
essentially use the Riemann-Lebesgue lemma, modified to take account of the dependence 
of the nonoscillatory terms on r. We have: 


r 9 _I2±2 r Q , TH. + 7 

I n = T 2 + 

where, following the usual transformation, 

= 2 Loo ~ G ( v + ]j^' T )) cos (^ nv ) dv - 

Here G is given by: 

G{v,t) = e' v2 (u m+3 d 1 (u,r) + u m+ V 2 (u, r)). 

Again G is absolutely integrable uniformly in r so we need only show that 

Bm(6>,T)-G(„ + ^, T )) = 0. 

Making use ol the fact that lor fixed v the limit r — > 00 implies u — ► 0 we have: 


(114) 


(115) 


(116) 


(117) 


Jim dj(u,r) = lim 




u — *-0 d* 


IT 


= -F 3 H 0 , 


and 


lim r) = lim (^Md ± g - jjgQ _ 

r — 00 u m+1 “ 1 


(118) 


(119) 


Fiom these we may conclude that (117) holds and, therefore, that lim T _ >00 7 3 = 0. We 
have shown, then, that I* = o(r-^). This, in turn, implies: 


/„ = o {e rJ ^T , t — ► 00 . 


(120) 

As the bounds obtained above may be made independent of n, we conclude that the 
contribution of the oscillatory terms is dominated by the contribution of the first term 
in the Fourier series, completing the proof of the theorem. 

It should be noted that we have not computed the leading order asymptotics of /„. 
One might at first glance expect that I\ determines the leading order behavior. We 
have, however, shown that it is exponentially small while the bounds obtained for the 
remaining terms decay algebraically. 
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